MODELOS LINEALES

Los modelos lineales son ampliamente utilizados en acuicultura para explicar, modelar o predecir la relaciĆ³n lineal de una variable respuesta \(Y\) con una o mĆ”s variables predictoras \(X_1,ā€¦,X_p\).

\[Y_{i} = \beta_{0} + \beta_{1} X_{i1} + \beta_{2} \beta_{i2} + ... + \beta_{p} X_{ip} + \epsilon_{i}\]

Supuestos: \(Y_i\) se distribuye normalmente

  • Si p = 1, el modelo es una regresiĆ³n lineal simple.

  • Si p > 1, el modelo es una regresiĆ³n lineal mĆŗltiple.

  • Si p > 1 y alguna variable predictora es categĆ³rica, el modelo se denomina ANCOVA.

PRUEBA DE HIPƓTESIS: COEFICIENTE DE REGRESIƓN E INTERCEPTO

Prueba de hipĆ³tesis del coeficiente de regresiĆ³n y el intercepto

La hipĆ³tesis nula en ambos casos es que el coeficiente del intercepto (\(Ī²_0\)) y de la pendiente (\(Ī²_1\)) son iguales a 0.

\(H_0:Ī²_0 = 0\) y \(H_0:Ī²_1 = 0\)

PREDICCIƓN A PARTIR DEL MODELO DE REGRESIƓN LINEAL SIMPLE

El Coeficiente de determinaciĆ³n (RĀ²) se calcula como el cuadrado del coeficiente de correlaciĆ³n de pearson. Este nos indica que tan buena es la predicciĆ³n que el modelo hace de los datos.

El R cuadrado ajustado (RĀ² ajustado) nos dice quĆ© porcentaje de la variaciĆ³n de la variable dependiente es explicado por la o las variables independientes de manera conjunta. En un modelo de regresiĆ³n lineal, el RĀ² ajustado es una medida de bondad que considera el nĆŗmero de variables existentes en el modelo.

\[R^2_{ajust} =1-(1-R^2)\frac{n-1}{n-p-1}\]

donde:

\(n\) = tamaƱo de la muestra

\(p\) = cantidad de variables predictoras en el modelo

Prueba de hipĆ³tesis del modelo completo

La hipĆ³tesis nula es si los coeficientes son iguales a 0.

\(H_0:Ī²_j = 0\) ; \(j = 1, 2,...,k\)

Objetivos de aprendizaje

Los objetivos de aprendizaje de esta guĆ­a son:

1. - Realizar anƔlisis de datos usando modelos lineales.

2. - Realizar grƔficas avanzadas con ggplot2.

3. - Elaborar un reporte dinƔmico en formato html con Rmarkdown.

EJERCICIOS

ESTUDIO DE CASO: MUTACIƓN GEN MIOSTATINA

En vertebrados, la miostatina (MSTN) actua como es un represor del crecimiento muscular. Muchos trabajos han descrito que mutacioones en ese gen producen un fenotipo llamado de doble mĆŗsculo caracterizado por un crecimiento exacerbado del mĆŗsculo esquelĆ©tico. En este ejercicio trabajaremos con un set de datos simulado de perros de carrera que poseen este distintivo fenotipo Adaptado de Mosher et al.Ā 2007

Figura 1: ComparaciĆ³n de animales con diferentes genotipos del gen Miostatina: (A) +/+; (B) 1 mutante con alelo cys ā†’ stop (mh/+); (C) Dos copias de la mutaciĆ³n cys ā†’ stop (mh/mh).

Figura 1: ComparaciĆ³n de animales con diferentes genotipos del gen Miostatina: (A) +/+; (B) 1 mutante con alelo cys ā†’ stop (mh/+); (C) Dos copias de la mutaciĆ³n cys ā†’ stop (mh/mh).

La variable respuesta se denomina mass-to-height ratio y se expresa en Kg/cm. Para este set de datos existen dos variables predictoras que se deberƔn analizar: 1) Geno: que corresponde al genotipo de los animales; 2) Sex: que corresponde al sexo.

La variable Geno tiene 3 niveles denominados 0, 1 y 2 equivalentes a los genotipos +/+, +/mh y mh/mh donde 0, 1 y 2 representan el nĆŗmero de alelos mh. Esta forma numĆ©rica de representar los genotipos es conveniente para analizar la variable Geno como una variable cuantitativa discreta 0, 1 y 2 o como un factor (ej. factor(Geno)) segĆŗn sea nuestro interĆ©s. La variable Sex tiene dos niveles male y female.

Ejercicio 1. Elaborar y configurar Reporte en formato .Rmd

Elabore un documento .Rmd y configure su reporte para exportar en .html. Instale los paquetes readxl, ggplot2, dplyr y multcomp para el anƔlisis de los datos.

knitr::opts_chunk$set(echo = TRUE)
library(readxl)
library(ggplot2)
library(dplyr)
library(multcomp)

Ejecute cada uno de los siguientes ejercicios en uno o mĆ”s bloques de cĆ³digos diferentes. Sea ordenado y documente su reporte adecuadamente.

Ejercicio 2. Exploratorio set de datos snp.data

Importe el set de datos myostatin.deletion.data.xlsx y transforme solo la variable Sex a factor, luego realice un anƔlisis exploratorio de datos.

Incluya:

a). Resumen estadĆ­stico de todas las variables.

snp.data <- read_excel("myostatin.deletion.data.xlsx")

snp.data$Sex <- as.factor(snp.data$Sex)
summary(snp.data)
##      animal         Geno            Sex      mass-to-height ratio
##  Min.   :  1   Min.   :0.0000   female:494   Min.   :0.2009      
##  1st Qu.:174   1st Qu.:0.0000   male  :199   1st Qu.:0.2919      
##  Median :347   Median :1.0000                Median :0.3401      
##  Mean   :347   Mean   :0.9769                Mean   :0.3590      
##  3rd Qu.:520   3rd Qu.:2.0000                3rd Qu.:0.4100      
##  Max.   :693   Max.   :2.0000                Max.   :0.5993

b). Tabla de frecuencia combinada Geno por Sex, esto es importante para ver si los datos son o no balanceados.

table(snp.data$Sex, snp.data$Geno)
##         
##            0   1   2
##   female 155  80 259
##   male   120  79   0
tapply(snp.data$`mass-to-height ratio`, factor(snp.data$Geno), mean)
##         0         1         2 
## 0.2829893 0.3263597 0.4597127
tapply(snp.data$`mass-to-height ratio`, factor(snp.data$Sex), mean)
##    female      male 
## 0.3756034 0.3177426

c). Histograma de la variable respuesta mass-to-height ratio.

ggplot(snp.data, aes(x=`mass-to-height ratio`))+
  geom_histogram(color="darkblue", fill="lightblue", bins = 12)

Ejercicio 3. RegresiĆ³n lineal simple.

a). Realice una grĆ”fica de dispersiĆ³n (geom_point) de weight en funciĆ³n de Geno, incluya el comando geom_smooth(method=lm) para agregar la lĆ­nea de regresiĆ³n a la grĆ”fica.

p <- ggplot(snp.data, aes(x = Geno, y = `mass-to-height ratio`))
p + geom_point() + xlab("Reference allele count") + geom_smooth(method=lm)

b). Realice un anĆ”lisis de regresiĆ³n lineal para investigar la asociaciĆ³n entre weight y Geno usando las funciones lm(), summary(). Interprete los resultados del modelo lineal y responda las siguientes preguntas. ĀæQuĆ© representa el estimador de Geno y del intercepto?. ĀæLa pendiente es significativamente distinta de cero?. ĀæY el intercepto?

lm.geno <- lm(`mass-to-height ratio` ~ Geno, data = snp.data)
summary(lm.geno)
## 
## Call:
## lm(formula = `mass-to-height ratio` ~ Geno, data = snp.data)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.110802 -0.044236 -0.006808  0.044296  0.150224 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 0.272969   0.003434   79.49   <2e-16 ***
## Geno        0.088052   0.002615   33.67   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.06041 on 691 degrees of freedom
## Multiple R-squared:  0.6213, Adjusted R-squared:  0.6208 
## F-statistic:  1134 on 1 and 691 DF,  p-value: < 2.2e-16

c). Realice la predicciĆ³n de la media poblacional de weight para los genotipos AA, AB y BB y la predicciĆ³n del peso weight de tres nuevos individuos con genotipo AA, AB y BB. ĀæQuĆ© tan buena son las predicciones obtenidas? Analice el RĀ² ajustado y los intervalos de confianza.

predict.lm(lm.geno, newdata=data.frame(Geno=c(0,1,2)), interval="confidence")
##         fit       lwr       upr
## 1 0.2729690 0.2662268 0.2797112
## 2 0.3610212 0.3565142 0.3655282
## 3 0.4490734 0.4421530 0.4559937
predict.lm(lm.geno, newdata=data.frame(Geno=c(0,1,2)), interval="prediction")
##         fit       lwr       upr
## 1 0.2729690 0.1541725 0.3917655
## 2 0.3610212 0.2423306 0.4797118
## 3 0.4490734 0.3302666 0.5678801

d). Realice una anova del modelo lineal con la funciĆ³n anova() (supone varianzas homogeneas) y compare con una anova que no supone varianzas homogeneas oneway.test()

# anova suponiendo varianzas iguales
anova(lm.geno)
## Analysis of Variance Table
## 
## Response: mass-to-height ratio
##            Df Sum Sq Mean Sq F value    Pr(>F)    
## Geno        1 4.1373  4.1373  1133.8 < 2.2e-16 ***
## Residuals 691 2.5215  0.0036                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# anova suponiendo varianzas diferentes
oneway.test(`mass-to-height ratio` ~ Geno, data = snp.data)
## 
##  One-way analysis of means (not assuming equal variances)
## 
## data:  `mass-to-height ratio` and Geno
## F = 520.26, num df = 2.00, denom df = 427.93, p-value < 2.2e-16

e) El siguiente modelo lineal permite realizar pruebas de contraste entre los tres diferentes genotipos.

lm(weight ~ -1 + factor(Geno), data = snp.data)

Compare los Betas y el \[R^2\] de este modelo con el realizado en el ejercicio b.

lm.geno2 <- lm(`mass-to-height ratio` ~ -1 + factor(Geno), data = snp.data)

summary(lm.geno2)
## 
## Call:
## lm(formula = `mass-to-height ratio` ~ -1 + factor(Geno), data = snp.data)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.119005 -0.044944 -0.000312  0.040813  0.139584 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## factor(Geno)0 0.282989   0.003461   81.75   <2e-16 ***
## factor(Geno)1 0.326360   0.004552   71.69   <2e-16 ***
## factor(Geno)2 0.459713   0.003567  128.89   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.0574 on 690 degrees of freedom
## Multiple R-squared:  0.9763, Adjusted R-squared:  0.9762 
## F-statistic:  9478 on 3 and 690 DF,  p-value: < 2.2e-16
summary(lm.geno)
## 
## Call:
## lm(formula = `mass-to-height ratio` ~ Geno, data = snp.data)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.110802 -0.044236 -0.006808  0.044296  0.150224 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 0.272969   0.003434   79.49   <2e-16 ***
## Geno        0.088052   0.002615   33.67   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.06041 on 691 degrees of freedom
## Multiple R-squared:  0.6213, Adjusted R-squared:  0.6208 
## F-statistic:  1134 on 1 and 691 DF,  p-value: < 2.2e-16

f) Investigue y aplique los comandos contrMat() y glht para realizar una prueba de contrastes con ajuste de bonferroni para comparaciones mĆŗltiples.

# Elabora matriz de contrastes para el factor Geno
contrastes = contrMat(table(snp.data$Geno), type="Tukey")

# Realiza comparaciones multiples
mc2 = glht(lm.geno2, linfct =contrastes)
summary(mc2, test=adjusted("bonferroni"))
## 
##   Simultaneous Tests for General Linear Hypotheses
## 
## Multiple Comparisons of Means: Tukey Contrasts
## 
## 
## Fit: lm(formula = `mass-to-height ratio` ~ -1 + factor(Geno), data = snp.data)
## 
## Linear Hypotheses:
##            Estimate Std. Error t value Pr(>|t|)    
## 1 - 0 == 0 0.043370   0.005719   7.584 3.26e-13 ***
## 2 - 0 == 0 0.176723   0.004970  35.556  < 2e-16 ***
## 2 - 1 == 0 0.133353   0.005783  23.059  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- bonferroni method)

Ejercicio 4. AnƔlisis de covarianza

a). Realice una grĆ”fica de dispersiĆ³n (geom_point) de mass-to-height ratio en funciĆ³n de Geno y Sex, incluya el comando geom_smooth(method=lm) para agregar la lĆ­nea de regresiĆ³n a la grĆ”fica.

q <- ggplot(snp.data, aes(x = Geno, y = `mass-to-height ratio`, shape=Sex, color=Sex))
q + geom_point() + xlab("Reference allele count") + geom_smooth(method=lm)

b). Realice un anĆ”lisis de covarianza para investigar la asociaciĆ³n entre mass-to-height ratio, Geno y Sex usando las funciones lm(), summary().

lm.geno.sex <- lm(`mass-to-height ratio` ~ Geno + Sex, data = snp.data)
summary(lm.geno.sex)
## 
## Call:
## lm(formula = `mass-to-height ratio` ~ Geno + Sex, data = snp.data)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.107264 -0.046971 -0.006527  0.041388  0.151325 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 0.264638   0.004396  60.205  < 2e-16 ***
## Geno        0.091667   0.002864  32.006  < 2e-16 ***
## Sexmale     0.016714   0.005555   3.009  0.00272 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.06006 on 690 degrees of freedom
## Multiple R-squared:  0.6262, Adjusted R-squared:  0.6251 
## F-statistic:   578 on 2 and 690 DF,  p-value: < 2.2e-16

c). Interprete los resultados del modelo lineal.

Ejercicio 5. Anova 2 vĆ­as con interacciĆ³n

a). Realice grĆ”ficas del tamaƱo de los efectos plot.design() y de interacciĆ³n interaction.plot de Geno y Sex sobre mass-to-height ratio.

plot.design(snp.data$`mass-to-height ratio` ~ factor(snp.data$Geno) + factor(snp.data$Sex))

interaction.plot(factor(snp.data$Geno), factor(snp.data$Sex), snp.data$`mass-to-height ratio`, mean)

b). Realice un anĆ”lisis de varianza con interacciĆ³n para investigar la asociaciĆ³n entre weight, Geno y Sex usando las funciones lm(), summary().

lm.qtl <- lm(`mass-to-height ratio` ~ factor(Geno) * factor(Sex), data = snp.data)
summary(lm.qtl)
## 
## Call:
## lm(formula = `mass-to-height ratio` ~ factor(Geno) * factor(Sex), 
##     data = snp.data)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.119005 -0.041228  0.000675  0.040398  0.139584 
## 
## Coefficients: (1 not defined because of singularities)
##                               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                   0.271157   0.004495  60.321  < 2e-16 ***
## factor(Geno)1                 0.034509   0.007704   4.479 8.78e-06 ***
## factor(Geno)2                 0.188556   0.005683  33.177  < 2e-16 ***
## factor(Sex)male               0.027117   0.006805   3.985 7.47e-05 ***
## factor(Geno)1:factor(Sex)male 0.014534   0.011185   1.299    0.194    
## factor(Geno)2:factor(Sex)male       NA         NA      NA       NA    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.05597 on 688 degrees of freedom
## Multiple R-squared:  0.6764, Adjusted R-squared:  0.6745 
## F-statistic: 359.5 on 4 and 688 DF,  p-value: < 2.2e-16

c). Interprete los resultados del modelo lineal.

Ejercicio 6. Evalue supuestos del modelo

a). Evalue mormalidad del modelo con interacciĆ³n del ejercicio 5 usando mĆ©todos grĆ”ficos y prueba de shapiro.

shapiro.test (residuals (lm.qtl))
## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(lm.qtl)
## W = 0.99161, p-value = 0.0005773
plot(lm.qtl, which = 2)

b). Evalue homogeneidad de varianzas usando el comando plot.

plot(lm.qtl, which = 1)